Research Question: Do damselfish (Dascyllus flavicaudus) accelerate coral wound healing rates in Pocillopora spp.?
Experimental Design:
This document contains:
# ==================== Core Packages ====================
library(tidyverse) # Data wrangling and visualization
library(here) # Relative file paths
library(janitor) # Data cleaning
# ==================== Statistical Modeling ====================
library(lme4) # Linear mixed-effects models
library(glmmTMB) # Generalized linear mixed models
library(DHARMa) # Residual diagnostics for GLMMs
library(MuMIn) # Model selection
library(parameters) # Model parameter extraction
library(broom.mixed) # Tidy model outputs
library(ggeffects) # Marginal effects
# ==================== Visualization ====================
library(ggpubr) # Publication-ready plots
library(viridis) # Color palettes
library(patchwork) # Combine plots
library(gt) # Publication-quality tables
# ==================== Fonts ====================
library(sysfonts) # Font management
library(showtext) # Font rendering
# ==================== Optional: Advanced Models ====================
if (!requireNamespace("gamlss", quietly = TRUE)) {
install.packages("gamlss")
}
library(gamlss)
library(gamlss.dist)
# Setup font rendering
font_add_google(name = "Josefin Sans", family = "josefin")
showtext_auto()
# Set default theme
theme_set(theme_bw() + theme(panel.grid = element_blank(), plot.background = element_blank()))# Import wound closure binary data (healed yes/no)
wound_closure <- read_csv(
here("data", "fish_regen_mastersheet_wound closure_necrosis_sa - mastersheet.csv")
) %>%
clean_names()
# Import quantitative healing rate data
percent_closure <- read_csv(
here("data", "final_wound_size_dascyllus_project - Sheet1.csv")
) %>%
clean_names()
cat("✓ Data imported successfully\n")## ✓ Data imported successfully
## - Wound closure binary data: 72 observations
## - Percent closure data: 48 observations
#' Convert numeric treatment codes to labeled factors
#'
#' @param df Data frame containing fish and wound columns
#' @return Data frame with properly ordered factor levels
prep_treatments <- function(df) {
df %>%
mutate(
# Convert fish coding: 1 = Fish present, 0 = Fish absent
fish = ifelse(fish == 1, "Fish", ifelse(fish == 0, "No Fish", fish)),
# Convert wound size coding: 0 = None, 1 = Small, 2 = Large
wound = case_when(
wound == 0 ~ "No Wound",
wound == 1 ~ "Small",
wound == 2 ~ "Large",
TRUE ~ as.character(wound)
)
) %>%
mutate(
# Set factor levels in logical order
fish = factor(fish, levels = c("No Fish", "Fish")),
wound = factor(wound, levels = c("No Wound", "Small", "Large"))
)
}# Process binary wound closure data
wound_closure_clean <- wound_closure %>%
filter(!is.na(wound_close), wound_close != "") %>%
# Standardize "barely healed" responses to "no"
mutate(wound_close = ifelse(wound_close %in% c("no-barely", "no - barely"), "no", wound_close)) %>%
mutate(wound_close = ifelse(wound_close == "yes", "Yes", "No")) %>%
prep_treatments() %>%
mutate(wound_close = as.character(wound_close))
# Create binary version for GLM modeling
wound_closure_glm <- wound_closure %>%
filter(!is.na(wound_close), wound_close != "") %>%
mutate(wound_bin = ifelse(wound_close == "yes", 1, 0)) %>%
prep_treatments() %>%
filter(!is.na(fish), !is.na(wound), !is.na(tank)) %>%
droplevels()
# Quality check: Ensure factors have adequate levels
if (nlevels(wound_closure_glm$fish) < 2 | nlevels(wound_closure_glm$wound) < 2) {
stop("ERROR: Fish or wound factors have <2 levels. Check input data.")
}
cat("✓ Data cleaning complete\n")## ✓ Data cleaning complete
## - Binary healing data: 48 corals
# Process continuous healing rate data
percent_closure <- percent_closure %>%
filter(percent_healed >= 0) %>% # Remove invalid percentages
prep_treatments() %>%
mutate(tank = as.factor(tank))
# Convert percent to proportion (0-1 scale)
percent_closure$prop <- percent_closure$percent_healed / 100
percent_closure$prop_remaining <- 1 - percent_closure$prop
cat("✓ Healing rate data prepared\n")## ✓ Healing rate data prepared
## - Valid observations: 47 corals
cat(" - Mean healing rate:", round(mean(percent_closure$healing_rate, na.rm = TRUE), 3), "cm²/day\n")## - Mean healing rate: 0.119 cm²/day
##
## === Treatment Sample Sizes ===
##
## No Wound Small Large
## No Fish 0 11 12
## Fish 0 12 12
# Calculate healing proportions by treatment
wound_prop <- wound_closure_clean %>%
group_by(wound, fish) %>%
summarize(
n = n(),
prop_healed = mean(wound_close == "Yes"),
.groups = "drop"
)
cat("\n=== Proportion Healed by Treatment ===\n")##
## === Proportion Healed by Treatment ===
## # A tibble: 4 × 4
## wound fish n prop_healed
## <fct> <fct> <int> <dbl>
## 1 Small No Fish 12 0.583
## 2 Small Fish 12 0.917
## 3 Large No Fish 12 0.417
## 4 Large Fish 12 0.75
# Summarize initial wound sizes
wounds_avgs <- percent_closure %>%
group_by(wound) %>%
summarize(
mean_size_cm2 = mean(initial_wound_size_cm, na.rm = TRUE),
sd_size_cm2 = sd(initial_wound_size_cm, na.rm = TRUE),
n = n(),
.groups = "drop"
) %>%
filter(wound != "No Wound")
# Display as formatted table
wounds_avgs %>%
gt() %>%
tab_header(title = "Initial Wound Size by Treatment") %>%
fmt_number(columns = c(mean_size_cm2, sd_size_cm2), decimals = 2) %>%
cols_label(
wound = "Wound Size",
mean_size_cm2 = "Mean (cm²)",
sd_size_cm2 = "SD (cm²)",
n = "N"
)| Initial Wound Size by Treatment | |||
| Wound Size | Mean (cm²) | SD (cm²) | N |
|---|---|---|---|
| Small | 1.48 | 0.40 | 23 |
| Large | 3.82 | 1.08 | 24 |
We fit linear mixed-effects models to test how fish presence and wound size affect healing rate (continuous response variable in cm²/day).
Model structure: - Fixed effects: Fish presence, Wound size, and their interaction - Random effect: Tank (accounts for shared tank environment) - Response: Healing rate (cm²/day)
# Full model with interaction
rate_full <- lmer(
healing_rate ~ fish * wound + (1 | tank),
data = percent_closure
)
# Additive model (no interaction)
rate_no_interaction <- lmer(
healing_rate ~ fish + wound + (1 | tank),
data = percent_closure
)
# Fish-only model
rate_no_fish <- lmer(
healing_rate ~ wound + (1 | tank),
data = percent_closure
)
# Wound-only model
rate_no_wound <- lmer(
healing_rate ~ fish + (1 | tank),
data = percent_closure
)
cat("✓ Models fitted successfully\n")## ✓ Models fitted successfully
# Test for interaction effect
lrt_interaction_zero_rate <- anova(rate_no_interaction, rate_full, test = "Chisq")
# Test for fish main effect
lrt_fish_zero_rate <- anova(rate_no_fish, rate_no_interaction, test = "Chisq")
# Test for wound main effect
lrt_wound_zero_rate <- anova(rate_no_wound, rate_no_interaction, test = "Chisq")
# Compile results into table
lrt_table_rate <- tibble(
Test = c("Interaction (Fish × Wound)", "Fish Effect", "Wound Effect"),
ChiSq = c(
lrt_interaction_zero_rate$Chisq[2],
lrt_fish_zero_rate$Chisq[2],
lrt_wound_zero_rate$Chisq[2]
),
Df = c(
lrt_interaction_zero_rate$Df[2],
lrt_fish_zero_rate$Df[2],
lrt_wound_zero_rate$Df[2]
),
p_value = c(
lrt_interaction_zero_rate$`Pr(>Chisq)`[2],
lrt_fish_zero_rate$`Pr(>Chisq)`[2],
lrt_wound_zero_rate$`Pr(>Chisq)`[2]
)
) %>%
mutate(across(where(is.numeric), round, 3))
# Format as publication table
gt_table_rate <- lrt_table_rate %>%
gt() %>%
tab_header(
title = "Likelihood Ratio Tests for Effects on Healing Rate"
) %>%
cols_label(
Test = "Model Comparison",
ChiSq = "χ²",
Df = "df",
p_value = "P-value"
) %>%
fmt_number(columns = c(ChiSq, p_value), decimals = 3) %>%
tab_options(
table.font.size = 14,
heading.title.font.size = 16
)
# Display table
print(gt_table_rate)## <div id="ldbicnuyvi" style="padding-left:0px;padding-right:0px;padding-top:10px;padding-bottom:10px;overflow-x:auto;overflow-y:auto;width:auto;height:auto;">
## <style>#ldbicnuyvi table {
## font-family: system-ui, 'Segoe UI', Roboto, Helvetica, Arial, sans-serif, 'Apple Color Emoji', 'Segoe UI Emoji', 'Segoe UI Symbol', 'Noto Color Emoji';
## -webkit-font-smoothing: antialiased;
## -moz-osx-font-smoothing: grayscale;
## }
##
## #ldbicnuyvi thead, #ldbicnuyvi tbody, #ldbicnuyvi tfoot, #ldbicnuyvi tr, #ldbicnuyvi td, #ldbicnuyvi th {
## border-style: none;
## }
##
## #ldbicnuyvi p {
## margin: 0;
## padding: 0;
## }
##
## #ldbicnuyvi .gt_table {
## display: table;
## border-collapse: collapse;
## line-height: normal;
## margin-left: auto;
## margin-right: auto;
## color: #333333;
## font-size: 14px;
## font-weight: normal;
## font-style: normal;
## background-color: #FFFFFF;
## width: auto;
## border-top-style: solid;
## border-top-width: 2px;
## border-top-color: #A8A8A8;
## border-right-style: none;
## border-right-width: 2px;
## border-right-color: #D3D3D3;
## border-bottom-style: solid;
## border-bottom-width: 2px;
## border-bottom-color: #A8A8A8;
## border-left-style: none;
## border-left-width: 2px;
## border-left-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_caption {
## padding-top: 4px;
## padding-bottom: 4px;
## }
##
## #ldbicnuyvi .gt_title {
## color: #333333;
## font-size: 16px;
## font-weight: initial;
## padding-top: 4px;
## padding-bottom: 4px;
## padding-left: 5px;
## padding-right: 5px;
## border-bottom-color: #FFFFFF;
## border-bottom-width: 0;
## }
##
## #ldbicnuyvi .gt_subtitle {
## color: #333333;
## font-size: 85%;
## font-weight: initial;
## padding-top: 3px;
## padding-bottom: 5px;
## padding-left: 5px;
## padding-right: 5px;
## border-top-color: #FFFFFF;
## border-top-width: 0;
## }
##
## #ldbicnuyvi .gt_heading {
## background-color: #FFFFFF;
## text-align: center;
## border-bottom-color: #FFFFFF;
## border-left-style: none;
## border-left-width: 1px;
## border-left-color: #D3D3D3;
## border-right-style: none;
## border-right-width: 1px;
## border-right-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_bottom_border {
## border-bottom-style: solid;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_col_headings {
## border-top-style: solid;
## border-top-width: 2px;
## border-top-color: #D3D3D3;
## border-bottom-style: solid;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## border-left-style: none;
## border-left-width: 1px;
## border-left-color: #D3D3D3;
## border-right-style: none;
## border-right-width: 1px;
## border-right-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_col_heading {
## color: #333333;
## background-color: #FFFFFF;
## font-size: 100%;
## font-weight: normal;
## text-transform: inherit;
## border-left-style: none;
## border-left-width: 1px;
## border-left-color: #D3D3D3;
## border-right-style: none;
## border-right-width: 1px;
## border-right-color: #D3D3D3;
## vertical-align: bottom;
## padding-top: 5px;
## padding-bottom: 6px;
## padding-left: 5px;
## padding-right: 5px;
## overflow-x: hidden;
## }
##
## #ldbicnuyvi .gt_column_spanner_outer {
## color: #333333;
## background-color: #FFFFFF;
## font-size: 100%;
## font-weight: normal;
## text-transform: inherit;
## padding-top: 0;
## padding-bottom: 0;
## padding-left: 4px;
## padding-right: 4px;
## }
##
## #ldbicnuyvi .gt_column_spanner_outer:first-child {
## padding-left: 0;
## }
##
## #ldbicnuyvi .gt_column_spanner_outer:last-child {
## padding-right: 0;
## }
##
## #ldbicnuyvi .gt_column_spanner {
## border-bottom-style: solid;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## vertical-align: bottom;
## padding-top: 5px;
## padding-bottom: 5px;
## overflow-x: hidden;
## display: inline-block;
## width: 100%;
## }
##
## #ldbicnuyvi .gt_spanner_row {
## border-bottom-style: hidden;
## }
##
## #ldbicnuyvi .gt_group_heading {
## padding-top: 8px;
## padding-bottom: 8px;
## padding-left: 5px;
## padding-right: 5px;
## color: #333333;
## background-color: #FFFFFF;
## font-size: 100%;
## font-weight: initial;
## text-transform: inherit;
## border-top-style: solid;
## border-top-width: 2px;
## border-top-color: #D3D3D3;
## border-bottom-style: solid;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## border-left-style: none;
## border-left-width: 1px;
## border-left-color: #D3D3D3;
## border-right-style: none;
## border-right-width: 1px;
## border-right-color: #D3D3D3;
## vertical-align: middle;
## text-align: left;
## }
##
## #ldbicnuyvi .gt_empty_group_heading {
## padding: 0.5px;
## color: #333333;
## background-color: #FFFFFF;
## font-size: 100%;
## font-weight: initial;
## border-top-style: solid;
## border-top-width: 2px;
## border-top-color: #D3D3D3;
## border-bottom-style: solid;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## vertical-align: middle;
## }
##
## #ldbicnuyvi .gt_from_md > :first-child {
## margin-top: 0;
## }
##
## #ldbicnuyvi .gt_from_md > :last-child {
## margin-bottom: 0;
## }
##
## #ldbicnuyvi .gt_row {
## padding-top: 8px;
## padding-bottom: 8px;
## padding-left: 5px;
## padding-right: 5px;
## margin: 10px;
## border-top-style: solid;
## border-top-width: 1px;
## border-top-color: #D3D3D3;
## border-left-style: none;
## border-left-width: 1px;
## border-left-color: #D3D3D3;
## border-right-style: none;
## border-right-width: 1px;
## border-right-color: #D3D3D3;
## vertical-align: middle;
## overflow-x: hidden;
## }
##
## #ldbicnuyvi .gt_stub {
## color: #333333;
## background-color: #FFFFFF;
## font-size: 100%;
## font-weight: initial;
## text-transform: inherit;
## border-right-style: solid;
## border-right-width: 2px;
## border-right-color: #D3D3D3;
## padding-left: 5px;
## padding-right: 5px;
## }
##
## #ldbicnuyvi .gt_stub_row_group {
## color: #333333;
## background-color: #FFFFFF;
## font-size: 100%;
## font-weight: initial;
## text-transform: inherit;
## border-right-style: solid;
## border-right-width: 2px;
## border-right-color: #D3D3D3;
## padding-left: 5px;
## padding-right: 5px;
## vertical-align: top;
## }
##
## #ldbicnuyvi .gt_row_group_first td {
## border-top-width: 2px;
## }
##
## #ldbicnuyvi .gt_row_group_first th {
## border-top-width: 2px;
## }
##
## #ldbicnuyvi .gt_summary_row {
## color: #333333;
## background-color: #FFFFFF;
## text-transform: inherit;
## padding-top: 8px;
## padding-bottom: 8px;
## padding-left: 5px;
## padding-right: 5px;
## }
##
## #ldbicnuyvi .gt_first_summary_row {
## border-top-style: solid;
## border-top-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_first_summary_row.thick {
## border-top-width: 2px;
## }
##
## #ldbicnuyvi .gt_last_summary_row {
## padding-top: 8px;
## padding-bottom: 8px;
## padding-left: 5px;
## padding-right: 5px;
## border-bottom-style: solid;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_grand_summary_row {
## color: #333333;
## background-color: #FFFFFF;
## text-transform: inherit;
## padding-top: 8px;
## padding-bottom: 8px;
## padding-left: 5px;
## padding-right: 5px;
## }
##
## #ldbicnuyvi .gt_first_grand_summary_row {
## padding-top: 8px;
## padding-bottom: 8px;
## padding-left: 5px;
## padding-right: 5px;
## border-top-style: double;
## border-top-width: 6px;
## border-top-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_last_grand_summary_row_top {
## padding-top: 8px;
## padding-bottom: 8px;
## padding-left: 5px;
## padding-right: 5px;
## border-bottom-style: double;
## border-bottom-width: 6px;
## border-bottom-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_striped {
## background-color: rgba(128, 128, 128, 0.05);
## }
##
## #ldbicnuyvi .gt_table_body {
## border-top-style: solid;
## border-top-width: 2px;
## border-top-color: #D3D3D3;
## border-bottom-style: solid;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_footnotes {
## color: #333333;
## background-color: #FFFFFF;
## border-bottom-style: none;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## border-left-style: none;
## border-left-width: 2px;
## border-left-color: #D3D3D3;
## border-right-style: none;
## border-right-width: 2px;
## border-right-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_footnote {
## margin: 0px;
## font-size: 90%;
## padding-top: 4px;
## padding-bottom: 4px;
## padding-left: 5px;
## padding-right: 5px;
## }
##
## #ldbicnuyvi .gt_sourcenotes {
## color: #333333;
## background-color: #FFFFFF;
## border-bottom-style: none;
## border-bottom-width: 2px;
## border-bottom-color: #D3D3D3;
## border-left-style: none;
## border-left-width: 2px;
## border-left-color: #D3D3D3;
## border-right-style: none;
## border-right-width: 2px;
## border-right-color: #D3D3D3;
## }
##
## #ldbicnuyvi .gt_sourcenote {
## font-size: 90%;
## padding-top: 4px;
## padding-bottom: 4px;
## padding-left: 5px;
## padding-right: 5px;
## }
##
## #ldbicnuyvi .gt_left {
## text-align: left;
## }
##
## #ldbicnuyvi .gt_center {
## text-align: center;
## }
##
## #ldbicnuyvi .gt_right {
## text-align: right;
## font-variant-numeric: tabular-nums;
## }
##
## #ldbicnuyvi .gt_font_normal {
## font-weight: normal;
## }
##
## #ldbicnuyvi .gt_font_bold {
## font-weight: bold;
## }
##
## #ldbicnuyvi .gt_font_italic {
## font-style: italic;
## }
##
## #ldbicnuyvi .gt_super {
## font-size: 65%;
## }
##
## #ldbicnuyvi .gt_footnote_marks {
## font-size: 75%;
## vertical-align: 0.4em;
## position: initial;
## }
##
## #ldbicnuyvi .gt_asterisk {
## font-size: 100%;
## vertical-align: 0;
## }
##
## #ldbicnuyvi .gt_indent_1 {
## text-indent: 5px;
## }
##
## #ldbicnuyvi .gt_indent_2 {
## text-indent: 10px;
## }
##
## #ldbicnuyvi .gt_indent_3 {
## text-indent: 15px;
## }
##
## #ldbicnuyvi .gt_indent_4 {
## text-indent: 20px;
## }
##
## #ldbicnuyvi .gt_indent_5 {
## text-indent: 25px;
## }
##
## #ldbicnuyvi .katex-display {
## display: inline-flex !important;
## margin-bottom: 0.75em !important;
## }
##
## #ldbicnuyvi div.Reactable > div.rt-table > div.rt-thead > div.rt-tr.rt-tr-group-header > div.rt-th-group:after {
## height: 0px !important;
## }
## </style>
## <table class="gt_table" data-quarto-disable-processing="false" data-quarto-bootstrap="false">
## <thead>
## <tr class="gt_heading">
## <td colspan="4" class="gt_heading gt_title gt_font_normal gt_bottom_border" style>Likelihood Ratio Tests for Effects on Healing Rate</td>
## </tr>
##
## <tr class="gt_col_headings">
## <th class="gt_col_heading gt_columns_bottom_border gt_left" rowspan="1" colspan="1" scope="col" id="Test">Model Comparison</th>
## <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1" scope="col" id="ChiSq">χ²</th>
## <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1" scope="col" id="Df">df</th>
## <th class="gt_col_heading gt_columns_bottom_border gt_right" rowspan="1" colspan="1" scope="col" id="p_value">P-value</th>
## </tr>
## </thead>
## <tbody class="gt_table_body">
## <tr><td headers="Test" class="gt_row gt_left">Interaction (Fish × Wound)</td>
## <td headers="ChiSq" class="gt_row gt_right">4.416</td>
## <td headers="Df" class="gt_row gt_right">1</td>
## <td headers="p_value" class="gt_row gt_right">0.036</td></tr>
## <tr><td headers="Test" class="gt_row gt_left">Fish Effect</td>
## <td headers="ChiSq" class="gt_row gt_right">8.047</td>
## <td headers="Df" class="gt_row gt_right">1</td>
## <td headers="p_value" class="gt_row gt_right">0.005</td></tr>
## <tr><td headers="Test" class="gt_row gt_left">Wound Effect</td>
## <td headers="ChiSq" class="gt_row gt_right">40.900</td>
## <td headers="Df" class="gt_row gt_right">1</td>
## <td headers="p_value" class="gt_row gt_right">0.000</td></tr>
## </tbody>
##
##
## </table>
## </div>
## Linear mixed model fit by REML ['lmerMod']
## Formula: healing_rate ~ fish * wound + (1 | tank)
## Data: percent_closure
##
## REML criterion at convergence: -140
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.72955 -0.53791 0.04203 0.56225 2.13584
##
## Random effects:
## Groups Name Variance Std.Dev.
## tank (Intercept) 6.678e-05 0.008172
## Residual 1.752e-03 0.041852
## Number of obs: 47, groups: tank, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.05682 0.01348 4.214
## fishFish 0.02166 0.01871 1.158
## woundLarge 0.07441 0.01748 4.257
## fishFish:woundLarge 0.05077 0.02444 2.077
##
## Correlation of Fixed Effects:
## (Intr) fshFsh wndLrg
## fishFish -0.721
## woundLarge -0.677 0.488
## fshFsh:wndL 0.484 -0.668 -0.715
##
## === Fixed Effects (Full Model) ===
## # A tibble: 4 × 7
## effect term estimate std.error statistic conf.low conf.high
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 fixed (Intercept) 0.0568 0.0135 4.21 0.0304 0.0832
## 2 fixed fishFish 0.0217 0.0187 1.16 -0.0150 0.0583
## 3 fixed woundLarge 0.0744 0.0175 4.26 0.0402 0.109
## 4 fixed fishFish:woundLarge 0.0508 0.0244 2.08 0.00287 0.0987
# Create publication-quality figure
healing_plot <- percent_closure %>%
ggplot(aes(x = wound, y = healing_rate, color = fish, shape = fish)) +
# Individual data points (jittered and dodged)
geom_jitter(
aes(color = fish, shape = fish),
position = position_jitterdodge(jitter.width = 0.1, dodge.width = 0.6),
size = 2.5,
alpha = 0.25,
shape = 16
) +
# Mean points
stat_summary(
fun = mean,
geom = "point",
size = 5,
position = position_dodge(width = 0.6)
) +
# Error bars (mean ± SD)
stat_summary(
fun.data = function(x) {
data.frame(
y = mean(x),
ymin = mean(x) - sd(x),
ymax = mean(x) + sd(x)
)
},
geom = "errorbar",
width = 0,
size = 1,
position = position_dodge(width = 0.6)
) +
# Color and shape scales
scale_color_manual(
values = fish_palette,
name = "Fish Presence",
labels = c("No Fish", "Fish")
) +
scale_shape_manual(
values = c(16, 18),
name = "Fish Presence",
labels = c("No Fish", "Fish")
)+
scale_y_continuous(limits = c(0, 0.35)) +
# Theme and aesthetics
theme_minimal(base_size = 16) +
theme(
panel.grid = element_blank(),
axis.title = element_text(size = 20, face = "bold"),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 16),
legend.position = "top",
legend.title = element_text(size = 18, face = "bold"),
legend.text = element_text(size = 16),
legend.key.size = unit(1.5, "lines"),
legend.margin = margin(0, 0, 10, 0),
legend.box.spacing = unit(0.5, "lines"),
axis.line.x = element_line(color = "black", size = 0.5),
axis.line.y = element_line(color = "black", size = 0.5)
) +
labs(
y = "Healing Rate (cm²/day)",
x = "Wound Size"
)
# Display plot
print(healing_plot)# Save as PDF (vector, print-quality)
ggsave(
filename = here("figures", "wound_closure", "summary_line_healing_rate_by_wound_fish.pdf"),
plot = healing_plot,
width = 7,
height = 5,
dpi = 600,
device = cairo_pdf
)
# Save as PNG (raster, screen-quality)
ggsave(
filename = here("figures", "wound_closure", "summary_line_healing_rate_by_wound_fish.png"),
plot = healing_plot,
width = 8,
height = 6,
dpi = 300,
bg = "white"
)
cat("✓ Figure saved to figures/wound_closure/\n")## ✓ Figure saved to figures/wound_closure/
# Wound size effect
rate_avg_wound <- percent_closure %>%
group_by(wound) %>%
summarize(
mean_rate = mean(healing_rate, na.rm = TRUE),
sd_rate = sd(healing_rate, na.rm = TRUE),
n = n(),
.groups = "drop"
)
cat("=== Healing Rate by Wound Size ===\n")## === Healing Rate by Wound Size ===
## # A tibble: 2 × 4
## wound mean_rate sd_rate n
## <fct> <dbl> <dbl> <int>
## 1 Small 0.0682 0.0199 23
## 2 Large 0.167 0.0669 24
# Fish presence effect
rate_avg_fish <- percent_closure %>%
group_by(fish) %>%
summarize(
mean_rate = mean(healing_rate, na.rm = TRUE),
sd_rate = sd(healing_rate, na.rm = TRUE),
n = n(),
.groups = "drop"
)
cat("\n=== Healing Rate by Fish Presence ===\n")##
## === Healing Rate by Fish Presence ===
## # A tibble: 2 × 4
## fish mean_rate sd_rate n
## <fct> <dbl> <dbl> <int>
## 1 No Fish 0.0957 0.0566 23
## 2 Fish 0.141 0.0759 24
# Interaction: Fish × Wound
rate_avg_interaction <- percent_closure %>%
group_by(fish, wound) %>%
summarize(
mean_rate = mean(healing_rate, na.rm = TRUE),
sd_rate = sd(healing_rate, na.rm = TRUE),
n = n(),
.groups = "drop"
)
cat("\n=== Healing Rate by Fish × Wound ===\n")##
## === Healing Rate by Fish × Wound ===
## # A tibble: 4 × 5
## fish wound mean_rate sd_rate n
## <fct> <fct> <dbl> <dbl> <int>
## 1 No Fish Small 0.0570 0.0114 11
## 2 No Fish Large 0.131 0.0586 12
## 3 Fish Small 0.0785 0.0208 12
## 4 Fish Large 0.204 0.0553 12
library(DHARMa)
# Simulate residuals for full model
sim_rate_full <- simulateResiduals(fittedModel = rate_full, n = 1000)
# Save DHARMa diagnostic plot
png(
here("figures", "diagnostics", "sim_rate_full.png"),
width = 7, height = 6, units = "in", res = 600, bg = "white"
)
plot(sim_rate_full, main = "DHARMa Diagnostics: Full Model")
dev.off()## quartz_off_screen
## 2
##
## === Dispersion Test ===
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.92535, p-value = 0.754
## alternative hypothesis: two.sided
##
## === Zero-Inflation Test ===
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = NaN, p-value = 1
## alternative hypothesis: two.sided
model_list <- list(
Full = rate_full,
NoInteraction = rate_no_interaction,
NoFish = rate_no_fish,
NoWound = rate_no_wound
)
for (nm in names(model_list)) {
cat("\n========== Diagnostics for", nm, "model ==========\n")
sim <- simulateResiduals(model_list[[nm]], n = 1000)
# Save DHARMa plot
png(
here("figures", "diagnostics", paste0("resid_", nm, ".png")),
width = 7, height = 6, units = "in", res = 600, bg = "white"
)
plot(sim, main = paste("DHARMa Diagnostics:", nm, "Model"))
dev.off()
# Print test results
cat("\n Dispersion test:\n")
print(testDispersion(sim))
cat("\n Outlier test:\n")
print(testOutliers(sim))
}##
## ========== Diagnostics for Full model ==========
##
## Dispersion test:
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.92535, p-value = 0.754
## alternative hypothesis: two.sided
##
##
## Outlier test:
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim
## outliers at both margin(s) = 0, observations = 47, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.00000000 0.07548573
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0
##
##
## ========== Diagnostics for NoInteraction model ==========
##
## Dispersion test:
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.94733, p-value = 0.842
## alternative hypothesis: two.sided
##
##
## Outlier test:
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim
## outliers at both margin(s) = 1, observations = 47, p-value = 0.08972
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0005385317 0.1129377171
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.0212766
##
##
## ========== Diagnostics for NoFish model ==========
##
## Dispersion test:
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.97263, p-value = 0.99
## alternative hypothesis: two.sided
##
##
## Outlier test:
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim
## outliers at both margin(s) = 1, observations = 47, p-value = 0.08972
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.0005385317 0.1129377171
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0.0212766
##
##
## ========== Diagnostics for NoWound model ==========
##
## Dispersion test:
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.97448, p-value = 0.934
## alternative hypothesis: two.sided
##
##
## Outlier test:
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sim
## outliers at both margin(s) = 0, observations = 47, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
## 0.00000000 0.07548573
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 )
## 0
##
## ✓ All diagnostic plots saved to figures/diagnostics/
# Add predictions to data
percent_closure$pred_full <- predict(rate_full)
# Create scatterplot
p_fit <- ggplot(percent_closure, aes(x = pred_full, y = healing_rate)) +
geom_point(color = "steelblue", size = 2, alpha = 0.7) +
geom_smooth(method = "lm", se = TRUE, color = "darkred") +
labs(
x = "Predicted Healing Rate",
y = "Observed Healing Rate",
title = "Model Fit Check: rate_full"
) +
theme_bw(base_size = 14)
print(p_fit)library(patchwork)
# Residuals vs index
p_resid_plot <- ggplot(
data.frame(resid = residuals(rate_full, type = "pearson")),
aes(x = seq_along(resid), y = resid)
) +
geom_point(alpha = 0.6) +
geom_hline(yintercept = 0, linetype = "dashed") +
theme_minimal() +
labs(title = "Pearson Residuals", x = "Observation", y = "Residual")
# Residuals vs fitted
p_fit_plot <- ggplot(
data.frame(fitted = fitted(rate_full),
resid = residuals(rate_full, type = "pearson")),
aes(x = fitted, y = resid)
) +
geom_point(alpha = 0.6) +
geom_hline(yintercept = 0, linetype = "dashed") +
theme_minimal() +
labs(title = "Residuals vs Fitted", x = "Fitted values", y = "Residuals")
# Combine panels
p_resid_panel <- (p_resid_plot | p_fit_plot) +
plot_annotation(
title = "Residual Diagnostics for Healing Rate Models",
theme = theme(plot.title = element_text(size = 18, face = "bold"))
)
print(p_resid_panel)p_qq <- ggplot(
data.frame(resid = residuals(rate_full, type = "pearson")),
aes(sample = resid)
) +
stat_qq(color = "steelblue", size = 2) +
stat_qq_line(color = "darkred", linewidth = 1) +
theme_bw(base_size = 14) +
labs(
title = "Normal Q-Q Plot of Residuals",
x = "Theoretical Quantiles",
y = "Sample Quantiles"
)
print(p_qq)p_scale_loc <- ggplot(
data.frame(
fitted = fitted(rate_full),
sqrt_abs_resid = sqrt(abs(residuals(rate_full, type = "pearson")))
),
aes(x = fitted, y = sqrt_abs_resid)
) +
geom_point(alpha = 0.6, color = "steelblue", size = 2) +
geom_smooth(se = TRUE, color = "darkred", method = "loess") +
theme_bw(base_size = 14) +
labs(
title = "Scale-Location Plot",
x = "Fitted Values",
y = expression(sqrt("|Standardized Residuals|"))
)
print(p_scale_loc)##
## === Model Diagnostics Summary ===
##
## Residual standard deviation: 0.04185224
## Number of observations: 47
##
## Residual summary statistics:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.114238 -0.022513 0.001759 0.000000 0.023531 0.089390
##
##
## === Checking for Influential Observations ===
influence_measures <- influence(rate_full, obs = TRUE)
cat("Cook's distance threshold (4/n):", 4/nobs(rate_full), "\n")## Cook's distance threshold (4/n): 0.08510638
## R version 4.5.1 (2025-06-13)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.0.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
##
## time zone: America/Los_Angeles
## tzcode source: internal
##
## attached base packages:
## [1] parallel splines stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] gamlss_5.5-0 nlme_3.1-168 gamlss.dist_6.1-1
## [4] gamlss.data_6.0-6 showtext_0.9-7 showtextdb_3.0
## [7] sysfonts_0.8.9 gt_1.0.0 patchwork_1.3.1
## [10] viridis_0.6.5 viridisLite_0.4.2 ggpubr_0.6.1
## [13] ggeffects_2.3.1 broom.mixed_0.2.9.6 parameters_0.28.0
## [16] MuMIn_1.48.11 DHARMa_0.4.7 glmmTMB_1.1.12
## [19] lme4_1.1-37 Matrix_1.7-3 janitor_2.2.1
## [22] here_1.0.1 lubridate_1.9.4 forcats_1.0.0
## [25] stringr_1.5.1 dplyr_1.1.4 purrr_1.1.0
## [28] readr_2.1.5 tidyr_1.3.1 tibble_3.3.0
## [31] ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] Rdpack_2.6.4 gridExtra_2.3 sandwich_3.1-1
## [4] rlang_1.1.6 magrittr_2.0.3 multcomp_1.4-28
## [7] snakecase_0.11.1 furrr_0.3.1 compiler_4.5.1
## [10] mgcv_1.9-3 systemfonts_1.2.3 vctrs_0.6.5
## [13] crayon_1.5.3 pkgconfig_2.0.3 fastmap_1.2.0
## [16] backports_1.5.0 labeling_0.4.3 utf8_1.2.6
## [19] rmarkdown_2.29 tzdb_0.5.0 nloptr_2.2.1
## [22] ragg_1.4.0 bit_4.6.0 xfun_0.53
## [25] cachem_1.1.0 jsonlite_2.0.0 broom_1.0.9
## [28] gap.datasets_0.0.6 R6_2.6.1 bslib_0.9.0
## [31] stringi_1.8.7 RColorBrewer_1.1-3 parallelly_1.45.1
## [34] car_3.1-3 boot_1.3-31 jquerylib_0.1.4
## [37] numDeriv_2016.8-1.1 estimability_1.5.1 Rcpp_1.1.0
## [40] knitr_1.50 zoo_1.8-14 timechange_0.3.0
## [43] tidyselect_1.2.1 abind_1.4-8 yaml_2.3.10
## [46] TMB_1.9.17 codetools_0.2-20 curl_7.0.0
## [49] listenv_0.9.1 lattice_0.22-7 withr_3.0.2
## [52] bayestestR_0.16.1 coda_0.19-4.1 evaluate_1.0.4
## [55] future_1.67.0 survival_3.8-3 xml2_1.4.0
## [58] pillar_1.11.0 gap_1.6 carData_3.0-5
## [61] stats4_4.5.1 reformulas_0.4.1 insight_1.4.0
## [64] generics_0.1.4 vroom_1.6.5 rprojroot_2.1.0
## [67] hms_1.1.3 scales_1.4.0 minqa_1.2.8
## [70] globals_0.18.0 xtable_1.8-4 glue_1.8.0
## [73] emmeans_1.11.2 tools_4.5.1 ggsignif_0.6.4
## [76] fs_1.6.6 mvtnorm_1.3-3 grid_4.5.1
## [79] rbibutils_2.3 datawizard_1.2.0 Formula_1.2-5
## [82] cli_3.6.5 textshaping_1.0.1 gtable_0.3.6
## [85] rstatix_0.7.2 sass_0.4.10 digest_0.6.37
## [88] TH.data_1.1-3 farver_2.1.2 htmltools_0.5.8.1
## [91] lifecycle_1.0.4 bit64_4.6.0-1 MASS_7.3-65
Analysis Complete ✓
All figures and tables saved to figures/ directory.